MATH 123 - Mathematical Modeling

Simulation Models

Instructor: Prof. Mario Bañuelos

# Import necessary packages
from IPython.display import HTML, Image
import numpy as np
# import plotting packages
import matplotlib.pyplot as plt
import plotly.graph_objects as go

HW Tips

If at time $t=0$, $I_0 = 5$ and then $I_1 = 15$. We can calculate the rate of infection by the following

$$ I_1 = I_0 + x S_0 I_0 $$

Solve for $x$ and use this as the new rate!

Outline

• Simulation Modeling
• Simulating Deterministic Behavior - Area Under Curve

Simulation Modeling

Up to now, we have been focusing on deterministic modeling. In case where we cannot explain a system analytically or data collected directly, we can simulate the behavior.

• Processes with an element of chance involved are called probabilistic, as opposed to deterministic, processes.

• Monte Carlo methods include a large class of computational algorithms that rely on repeated random sampling to obtain numerical results.

• Thus, these methods are probabilistic.

Simulating Deterministic Behavior - Area

Let us consider the equation of a circle

$$ x^2 + y^2 = 1 $$

and we wanted to calculate the area of the circle (without knowing the formula).

def area_circle(x,y,n):
    
    # initiate the count for inside the circle
    in_c = 0
    
    # check if the point (x,y) is inside the circle
    for i in range(n):
        if x[i]**2 + y[i]**2 <= 1:
            # increment the counter by 1
            # increase the number of points in the circle by 1
            in_c = in_c + 1
            #print('hey, (x,y) is in the circle')
    
    # calculate the fraction of points inside the circle
    in_cfrac = in_c / n
    
    # calculate the area of the quarter circle
    areac = in_cfrac
    
    return areac

$$ \frac{\text{area of circle}}{\text{area of square}} = \frac{\text{points inside the circle}}{\text{total points}} $$

# determine number of points
n = 10_000_000
# simulate x,y points (each point is inside [-1,1])
x = np.random.uniform(low=0, high=1, size=n)
y = np.random.uniform(low=0, high=1, size=n)

x.shape

(10000000,)

area_circle(x,y,n)*4

3.14175

# create/allocate memory for area 
my_area = np.zeros(4)
# fill each element with the calculated area
my_area[0] = area_circle(x,y,5)
my_area[1] = area_circle(x,y,10)
my_area[2] = area_circle(x,y,100)
my_area[3] = area_circle(x,y,1000)

fig = go.Figure()
fig.add_trace(go.Scatter(x=[5,10,100,1000], y = my_area, mode='markers'))

Simulating Deterministic Behavior (cont.)

We can do this for more than circles! If we wish to calculate the area under the curve of $y = f(x)$ over $ a \le x \le b$ by constructing a bounding rectangle with height $M$ and base length $b-a$.

Then,

$$ \frac{\text{area under curve}}{\text{area of rectangle}} \approx \frac{\text{number of points counted below curve}}{\text{total number of random points}} $$

Monte Carlo Area Algorithm

Input: Total number $n$ of random points to be generated in the simulation.

Output: Approximate area under $y=f(x)$ over $a \le x \le b$

1) Initialize: COUNTER = 0.

2) For $i = 1,2,\ldots,n$, do Steps 3–5.

3) Calculate random coordinates $x_i$ and $y_i$ that satisfy $a \le x_i \le b$ and $0\le y_i \le M$.

4) Calculate $f(x_i)$ for the random $x_i$ coordinate.

5) If $y_i \le f(x_i)$, then increment the COUNTER by 1. Otherwise, leave COUNTER as is.

6) Calculate AREA = $M(b - a)$ COUNTER/ n.

7) OUTPUT (AREA)

Practice

Use the Monte Carlo Area Algorithm to approximate the area under the curve of $y = \cos(x)$ from $-\pi/2 \le x \le \pi /2$.

np.pi

3.141592653589793

# np.cos() -- cosine
# np.pi -- pi

fig = go.Figure()
fig.add_trace(go.Scatter(x = x, y = np.cos(x), mode='markers'))
fig.show()

def area_cos(x,y,n):
    
    # initiate the count for inside the circle
    in_c = 0
    
    # check if the point (x,y) is below f(x)
    for i in range(n):
        if y[i] < np.cos(x[i]):
            # increment the counter by 1
            in_c = in_c + 1    
    # calculate the fraction of points below f(x)
    in_cfrac = in_c / n
    
    # calculate the area
    areac = 1 * np.pi * in_cfrac #[M]*(b-a)*in_cfrac
    
    return areac

# determine number of points
n = 10
# simulate x,y points (each point is inside [-1,1])
x = np.random.uniform(low=-np.pi/2, high=np.pi/2, size=n)
y = np.random.uniform(low=0, high=1, size=n)

area_cos(x,y,n)

2.5132741228718345

Practice

Use the Monte Carlo Area Algorithm to approximate the area under the curve of $y = x^3 - \sin(x)$.

def area_sin(x,y,n):
    
    # initiate the count for inside the circle
    in_c = 0
    
    # check if the point (x,y) is below f(x)
    for i in range(n):
        if y[i] < (x[i]**3) - np.sin(x[i]):
            # increment the counter by 1
            in_c = in_c + 1    
    # calculate the fraction of points below f(x)
    in_cfrac = in_c / n
    
    # calculate the area
    areac_new = 1 * np.pi * in_cfrac #[M]*(b-a)*in_cfrac
    
    return areac_new

n = 10
# simulate x,y points (each point is inside [-1,1])
x = np.random.uniform(low=-0.5, high=0, size=n)
y = np.random.uniform(low=0, high=1, size=n)

area_sin(x,y,n)

0.9424777960769379